// topologyTree.h: Build R-tree structure for fast search.
//
//////////////////////////////////////////////////////////////////////

#include "rtree/src/amr3bind.h"
#include "topologyTree.h"

//
//////////////////////////////////////////////////////////////////////////
topologyTree::topologyTree(void)
{
  m_nodeTree = NULL;
  m_elemTree = NULL;
  m_edgeTree = NULL;

  m_nodeIds = NULL;
  m_elemIds = NULL;
  m_edgeIds = NULL;

  m_mesh = NULL;

}

//
//////////////////////////////////////////////////////////////////////////
topologyTree::~topologyTree(void)
{
  if (m_nodeTree != NULL)
    RtreeDestroy (m_nodeTree);
  if (m_elemTree != NULL)
    RtreeDestroy (m_elemTree);
  if (m_edgeTree != NULL)
    RtreeDestroy (m_edgeTree);

  if (m_nodeIds != NULL)
    delete []m_nodeIds;
  if (m_elemIds != NULL)
    delete []m_elemIds;
  if (m_edgeIds != NULL)
    delete []m_edgeIds;
}

//
//////////////////////////////////////////////////////////////////////////
void topologyTree::setMesh( cTopology *mesh )
{
  m_mesh = mesh;
}

//
//////////////////////////////////////////////////////////////////////////
void topologyTree::buildNodeTree( double box )
{
  int i;
  double coord[3];

  // release current tree
  if (m_nodeTree != NULL)
    RtreeDestroy (m_nodeTree);
  if (m_nodeIds != NULL)
    delete []m_nodeIds;

  // create tree
  int nnodes = m_mesh->NumNodes();
  m_nodeIds = new int[nnodes];
  m_nodeTree = RtreeCreate ( );

  // build tree
  for (i = 0; i < nnodes; ++i)
  {
    m_nodeIds[i] = i;
    m_mesh->GetCoordNode(i, coord);
    xmin = coord[0] - box*0.5;
    ymin = coord[1] - box*0.5;
    zmin = coord[2] - box*0.5;
    xmax = coord[0] + box*0.5;
    ymax = coord[1] + box*0.5;
    zmax = coord[2] + box*0.5;
    RtreeInsert (m_nodeTree, (void*)&(m_nodeIds[i]), xmin, xmax, ymin, ymax, zmin, zmax);
  }
}

//
//////////////////////////////////////////////////////////////////////////
void topologyTree::buildElemTree()
{
  int i, k, nnodes, *conn;
  double coord[3];

  // release current tree
  if (m_elemTree != NULL)
    RtreeDestroy (m_elemTree);
  if (m_elemIds != NULL)
    delete []m_elemIds;

  // create tree
  int nelems = m_mesh->NumNodes();
  m_elemIds = new int[nelems];
  m_elemTree = RtreeCreate ( );


  // build tree
  for (i = 0; i < nelems; ++i)
  {
    nnodes = m_mesh->GetElemNNodes (i);
    if (nnodes == 0)
      continue;

    // init max and min
    conn = m_mesh->GetElemConn(i);
    m_mesh->GetCoordNode(conn[0], coord);
    xmin = xmax = coord[0];
    ymin = ymax = coord[1];
    zmin = zmax = coord[2];

    // get max and min
    for (k = 1; k < nnodes; ++k)
    {
      m_mesh->GetCoordNode(conn[k], coord);
      if (coord[0] < xmin)
        xmin = coord[0];
      if (coord[1] < ymin)
        ymin = coord[1];
      if (coord[2] < zmin)
        zmin = coord[2];
      if (coord[0] > xmax)
        xmax = coord[0];
      if (coord[1] > ymax)
        ymax = coord[1];
      if (coord[2] > zmax)
        zmax = coord[2];
    }

    // insert to elemTree
    m_elemIds[i] = i;
    RtreeInsert (m_elemTree, (void*)&(m_elemIds[i]), xmin, xmax, ymin, ymax, zmin, zmax);
  }

}

//
//////////////////////////////////////////////////////////////////////////
void topologyTree::buildEdgeTree()
{
  int i, ie, je;
  double coord[3];

  // release current tree
  if (m_edgeTree != NULL)
    RtreeDestroy (m_edgeTree);
  if (m_edgeIds != NULL)
    delete []m_edgeIds;

  // create tree
  int nedges = m_mesh->NumEdges();
  m_edgeIds = new int[nedges];
  m_edgeTree = RtreeCreate ( );


  // build tree
  for (i = 0; i < nedges; ++i)
  {
    m_mesh->GetEdge(i, &ie, &je);

    // init max and min
    m_mesh->GetCoordNode(ie, coord);
    xmin = xmax = coord[0];
    ymin = ymax = coord[1];
    zmin = zmax = coord[2];

    // get max and min
    m_mesh->GetCoordNode(je, coord);
    if (coord[0] < xmin)
      xmin = coord[0];
    if (coord[1] < ymin)
      ymin = coord[1];
    if (coord[2] < zmin)
      zmin = coord[2];
    if (coord[0] > xmax)
      xmax = coord[0];
    if (coord[1] > ymax)
      ymax = coord[1];
    if (coord[2] > zmax)
      zmax = coord[2];

    // insert to elemTree
    m_edgeIds[i] = i;
    RtreeInsert (m_edgeTree, (void*)&(m_edgeIds[i]), xmin, xmax, ymin, ymax, zmin, zmax);
  }
}

//
//////////////////////////////////////////////////////////////////////////
vector <int> topologyTree::findNodes
( double xmin, double ymin, double zmin, double xmax, double ymax, double zmax )
{
  return find (m_nodeTree, xmin, ymin, zmin, xmax, ymax, zmax);
}

//
//////////////////////////////////////////////////////////////////////////
vector <int> topologyTree::findElems
( double xmin, double ymin, double zmin, double xmax, double ymax, double zmax )
{
  return find (m_elemTree, xmin, ymin, zmin, xmax, ymax, zmax);
}

//
//////////////////////////////////////////////////////////////////////////
vector <int> topologyTree::findEdges
( double xmin, double ymin, double zmin, double xmax, double ymax, double zmax )
{
  return find (m_edgeTree, xmin, ymin, zmin, xmax, ymax, zmax);
}

//
//////////////////////////////////////////////////////////////////////////
vector <int> topologyTree::find( void *tree, double xmin, double ymin, double zmin, 
                                 double xmax, double ymax, double zmax)
{
  vector <int> res;
  void *inf;
  int   value;
  double xmn, xmx, ymn, ymx, zmn, zmx;

  if (tree != NULL)
  {
    RtreeInitSearchBox (tree, xmin, xmax, ymin, ymax, zmin, zmax);

    // get information
    while ((inf=RtreeSearchBox (tree, &xmn, &xmx, &ymn, &ymx, &zmn, &zmx)) != NULL)
    {
      value = *((int *) inf);
      res.push_back(value);
    }

  }

  return res;
}